R Code for GCB-19-0810

rm(list = ls())  # Remove all objects from memory

FIGURE 1. Illustration of the influence of most intense ENSO events on the coral cover in the ETP.

# Load packages 
library(dplyr)
library(RColorBrewer)
library(plotrix)
library(readr)
library(reshape)
library(data.table)

FIGURE 1a. ENSO and Coral Cover

# Plot ENOS 1971 - 2014
    ENOS<-read.table(file="http://www.cpc.ncep.noaa.gov/data/indices/ersst5.nino.mth.81-10.ascii",header=TRUE)
    ENOS1 <-ENOS[241:780, ]   # subset from 1982-2014
    Year <- ENOS1$YR
    Anom3 <- ENOS1$ANOM.3
    YR_ANOM3 <- data.frame(cbind(Year,Anom3))  # Create a data frame with two columns
    YR_ANOM3$Year <- as.numeric(YR_ANOM3$Year)
    x <- YR_ANOM3                              # Funcion My Plot to graph ENSO
    Year <- seq(
       as.Date(sprintf("%s/1/1",min(x$Year))), 
       as.Date(sprintf("%s/12/1",max(x$Year))), 
       "months")
   
     myPlot = function(y, x, 
                      colors=c("red","blue"), 
                      xlab=deparse(substitute(x)),
                      ylab=deparse(substitute(y)), ...)
    {
       n = length(x)
       y.pos = y.neg = y
       y.pos[y<0] = 0
       y.neg[y>0] = 0
       plot(y~x, type="n", xlab=xlab, ylab=ylab, ...)   # xaxt="n", Apague eje x
       #grid()
       polygon(c(x[1],x,x[n]),c(0,y.pos,0),col=colors[1],border=NA)
       polygon(c(x[1],x,x[n]),c(0,y.neg,0),col=colors[2],border=NA)
    } 

# Boxplot coral cover all years 
    cover <- read.csv("Data/Coral_Cover_ETP.csv", header=TRUE)      # Select file: Coral_Cover_ETP.csv
    cover$Coral_cover <- as.numeric(cover$Coral_cover)  # convert to numeric variable
    cover$Region <- as.factor(cover$Region)                # convert to nominal factors
    cover$Location <- as.factor(cover$Location)
    cover$Site <- as.factor(cover$Site) 
    cover$Year <- as.numeric(cover$Year)
    cover$Country <- factor(cover$Country) 

# Hierarchically aggregate by Site, Location and Region.
  # Agregate by site
    aggr.site <- aggregate(Coral_cover ~ Year+Region+Location+Site+Period_ENOS+Thermal_regime,FUN=mean,data=cover)
  # Agregate by location
    aggr.location <- aggregate(Coral_cover~Year+Region+Location+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.site)
  # Aggregate by region
    aggr.region <- aggregate(Coral_cover~Year+Region+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.location)

# Plot 
  # dev.off() # reset parameters for plot
  op <- par(mfcol=c(2,1), cex.main = 1.5, mar = c(0,5,0,2), oma=c(6,4,4,2), cex.lab = 1.5, cex.axis = 1,  las = 1)  # bty = "n",
  myPlot(x$A,Year,ylab="",xlab="",axes = FALSE)  
  axis(side=2, at=c(-2,-1,0,1,2), lwd = 0.5)
  box(lwd = 0.5)
  abline(h=0, lw = 0.4)  
  boxplot (aggr.region$Coral_cover ~ aggr.region$Year,  
           ylim = c(0, 100),           
           axes = FALSE,
           cex = 0.5,                 
           lwd = 1,
           pars = list(boxwex = 0.7, staplewex = 0.3, outwex = 0.3),
           las = 2,                  
           cex.axis=0.5,
           xlab = "Year",
           ylab = "") 
  axis(2,lwd = 0.5)
  box(lwd = 0.5)
  x <-c(1:45)
  YearsR = as.character(c(1970:2014))
  axis(1,  at = x, labels = YearsR, lwd = 0.5, cex.axis=0.5, las = 2)

  dev.off()
## null device 
##           1
  • Polynomial adjustment for all regions
library(ggplot2)

theme_set (theme_bw() + theme(panel.border = element_blank(), 
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(colour = "black")))
x = aggr.region$Year
y = aggr.region$Coral_cover
qplot(x,y, geom='smooth', span = 0.3)+ xlim(1970, 2014)+xlab("Year")+
   ylab("Coral cover")+ylim(0,80)+ geom_point()
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).

DHW and coral cover change

# Load packages 
library(tidyverse)
library(dplyr)
library(data.table)
library(RColorBrewer)

FIGURE 2A. The sequence of yearly maximum DHW experienced by the best-represented coral reefs sites from 1982 to 2014

  • Select and transform DHW data
# Select DHW data file: maxDHW.csv
  dhw <- read.csv("Data/maxDHW.csv", header=TRUE)

# Exclude years 1981, and 2017-2019
  dhw <- dhw[!dhw$Year == 1981, ] 
  dhw <- filter(dhw, Year < 2015)

# Create groups for plot
  ts.Cano_island <- subset(dhw, dhw$Site == "Cano_island")     
  ts.Uva_island <- subset(dhw, dhw$Site == "Uva_Island")          
  ts.Gorgona_island <- subset(dhw, dhw$Site == "Gorgona_Island")          
  ts.Darwin_island <- subset(dhw, dhw$Site == "Darwin_North_Anchorage")          
  ts.Bahia_Banderas <- subset(dhw, dhw$Site == "Bahia_Banderas")          
  • Plot max DHW in each location
# Year, maxDHW
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1)  
  par(lwd = 1)   
  
  plot(ts.Cano_island$Year,ts.Cano_island$maxDHW,   # Create a white false line with real data
       xlab="Years", ylab="Degree heating weeks (°C-weeks)",  cex=1.5, cex.lab=0.6, col="white",axes = FALSE,
       ylim = c(0, 28), lwd = 1.5) 
  par(new=T)
  rect(1982, 0, 1983.5, 28, col="#fee8c8", border ="#fee8c8")           # Strong ENSO years
  par(new=T)
  rect(1997, 0, 1998.5, 28, col="#fee8c8", border ="#fee8c8")
  par(new=T)
  plot(ts.Cano_island$Year, ts.Cano_island$maxDHW,type="l",              # Thermally stable
  xlab="", ylab="", lwd = 1.5, col="#8c2d04",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Uva_island$Year, ts.Uva_island$maxDHW, type="l",                # Thermally stable
       xlab="", ylab="", lwd = 1.5, col="#feb24c",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Gorgona_island$Year,ts.Gorgona_island$maxDHW, type="l",         # Tropical upwelling
       xlab="", ylab="",lwd = 1.5, col="#fb6a4a", axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Darwin_island$Year,ts.Darwin_island$maxDHW, type="l",           # Galapagos
       xlab="", ylab="", lwd = 1.5, col="#9e9ac8",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Bahia_Banderas$Year, ts.Bahia_Banderas$maxDHW, type="l",        # Seasonal
       xlab="", ylab="",lwd = 1.5,col="#a6d96a", axes = FALSE,
       ylim = c(0, 28))
  y <- c(0,4,8,12,16,20,24,28)                                            # Axis
  axis(2,lwd = 0.5, at = y,cex.axis=0.6)
  box(lwd = 0.5)
  x <-c(1982:2014)
  YearsR = as.character(c(1982:2014))
  axis(1,  at = x, labels = YearsR, lwd = 0.4, cex.axis=0.5, las = 2) 
  abline(h=4, col="#525252", lty=2)
  abline(h=8, col="#525252", lty=2)

# Legend
  sites <- c("Thermally stable - Caño Island","Thermally stable - Uva Island",
             "Tropical upwelling - Gorgona Island","Equatorial Upwelling - Darwin Island", 
             "Seasonal - Bahía Banderas")
  par(xpd=TRUE)     
  legend(x = 2000, y = 20, 
         xpd = NA,  # or TRUE 
         inset=c(-0.5,0),
         cex = 0.4,
         legend = sites,         # box categories 
         col = c("#8c2d04","#feb24c","#fb6a4a","#9e9ac8","#a6d96a"),
         lty= 1,
         lwd = 1.5,
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

  dev.off() # uncomment to reset parameters for plot 
## null device 
##           1

Relationship of coral cover rate of change and thermal stress

  • The DHW index measures the accumulated thermal stress above a bleaching threshold for 12 consecutive weeks.
  • At each site, the maximum monthly mean sea surface # temperature (MMM) is calculated from a long-term satellite climatology dataset.
  • Any temperature surpassing the MMM by 1ºC or more is added for 12 consecutive weeks,
  • because temperatures ~1ºC above the MMM are known to cause coral bleaching (Liu et al., 2006).
  • For example, if the water temperature during four consecutive weeks exceeds the MMM in 1.3°C, # 1.1°C, 1.2°C, and 1.4°C, then,
  • the site accumulates 5°C-weeks when these anomalies are summed.

We tested the hypothesis that the strength of the Log_annual_rate_of_change (response variable) is a function of the maxDHW experienced. We had measured the rate of change from four Thermal regimes, and fitted Thermal regimes as a random intercept, and estimate a common slope (change in the rate of change) for maxDHW across all Thermal regimes by fitting it as a fixed effect.

## Load data
  roc.cover <- read.csv("Data/RoC_DHW.csv", header=TRUE)  ### Select file:  RoC_DWH.csv

# Load packages
  library(lme4)
  library(lmerTest)
  library(MuMIn)
  • Model 6: coral cover change and DHW
model6 <- lmerTest::lmer(Log_annual_rate_of_change ~ maxDHW + (1|Thermal_regime), data=roc.cover)
step(model6)
## Backward reduced random-effect table:
## 
##                      Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>                             4 -17.738 43.476                     
## (1 | Thermal_regime)          0    3 -27.772 61.545 20.069  1  7.469e-06
##                         
## <none>                  
## (1 | Thermal_regime) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##        Eliminated  Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## maxDHW          0 0.65054 0.65054     1 128.59   9.873 0.002082 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
summary(model6) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
##    Data: roc.cover
## 
## REML criterion at convergence: 35.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7962 -0.1472  0.1167  0.4152  1.8621 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Thermal_regime (Intercept) 0.02501  0.1581  
##  Residual                   0.06589  0.2567  
## Number of obs: 133, groups:  Thermal_regime, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -0.085763   0.084256   3.155381  -1.018  0.38032   
## maxDHW       -0.016316   0.005193 128.594887  -3.142  0.00208 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## maxDHW -0.165
r.squaredGLMM(model6)   #  returns the marginal and the conditional R2
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.
##             R2m       R2c
## [1,] 0.05249512 0.3131505
# marginal R2  describes the proportion of variance explained by the fixed factor(s) alone:
# conditional R2 describes the proportion of variance explained by both the fixed and random factors
# Nakagawa, S. & Schielzeth, H. (2013). 

FIGURE 2B: Relationship of coral cover rate of change and thermal stress

# Load packages
  library(tidyverse)
  library(dplyr)
  library(data.table)
  library(RColorBrewer)
 
# Create groups for plot
  Seasonal          = subset(roc.cover, roc.cover$Thermal_regime == "Seasonal")          
  Thermally_Stable  = subset(roc.cover, roc.cover$Thermal_regime == "Thermally_Stable")  
  Tropical_Upwelling = subset(roc.cover, roc.cover$Thermal_regime == "Upwelling") 
  Galapagos         = subset(roc.cover, roc.cover$Thermal_regime == "Galapagos_Islands") 

# Plot
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1) 
  par(lwd = 1)
  plot (Seasonal$maxDHW, Seasonal$Log_annual_rate_of_change,
        pch = 21, lwd = 0.5, bg = "#a6d96a93", col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1),
        xlab = "Degree heating weeks (°C-weeks)",
        ylab = "Annual rate of change in coral cover")
  par(new=T); 
  plot (Thermally_Stable$maxDHW, Thermally_Stable$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, col="black", bg = "#feb24c93", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Tropical_Upwelling$maxDHW, Tropical_Upwelling$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#fb6a4a93",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Galapagos$maxDHW, Galapagos$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#9e9ac893",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  box(lwd = 0.5)
  x <-c(0,5,10,15,20,25,30)
  axis(1, lwd = 0.5, at = x, cex.axis = 0.5)
  y <- c(-1.5,-1.0,-0.5,0,0.5,1)
  axis(2,lwd = 0.5, at = y, cex.axis = 0.5)
  abline(v=4, col="#525252", lty=2); abline(v=8, col="#525252", lty=2); abline(h=0, col="#525252", lty=2); 
  # Axis text I, II, III, IV
  text(3.5, 0.3, "II", cex = 0.6); text(5, 0.3, "I", cex = 0.6)
  text(3.5, - 0.3, "III", cex = 0.6); text(5, -0.3, "IV", cex = 0.6)
  
  par(xpd=TRUE) 
  legend(x = 10, y = -0.5, 
         xpd = NA,   
         inset=c(-0.5,0),
         legend = c("Equatorial Upwelling", "Thermally stable", "Tropical upwelling","Seasonal"),
         pch=21,
         pt.bg = c("#9e9ac8","#fb6a4a","#fb6a4a","#a6d96a"),
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

  dev.off() # uncomment to reset parameters for plot
## null device 
##           1

Descriptives statistics heat stress

No stress

no.stress <- subset(roc.cover, maxDHW < 4)                      #subset of DHW < 4
length(no.stress$maxDHW)/length(roc.cover$maxDHW)               # ratio no stress: total              
## [1] 0.8496241
1- length(no.stress$maxDHW)/length(roc.cover$maxDHW)               # ratio stress: total              
## [1] 0.1503759
positive.rate <- subset(no.stress, Log_annual_rate_of_change > 0)   #subset of Rate > 0
length(positive.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 48
negative.rate <- subset(no.stress, Log_annual_rate_of_change < 0)   #subset of Rate > 0
length(negative.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 58

No stress seasonal

no.stress.seasonal <- subset(no.stress, Thermal_regime == "Seasonal")   
positive.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)   
## [1] 0.07142857

No stress upwelling

no.stress.upwelling <- subset(no.stress, Thermal_regime == "Upwelling")
positive.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

No stress thermally stable

no.stress.Thermally_Stable <- subset(no.stress, Thermal_regime == "Thermally_Stable")   
positive.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1.2

No stress Galapagos

no.stress.Galapagos_Islands <- subset(no.stress, Thermal_regime == "Galapagos_Islands") 
positive.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

Supplementary figures

library(pals)
library(dplyr)
library(reshape2)
library(RColorBrewer) 

FIGURE S2: a) Number of regions surveyed per year, b): Number of countries surveyed per year

# dev.off() # uncomment to reset parameters for plot

# Barplot number of surveys per year
  obs_num = as.data.frame(table(cover$Year,cover$Region))    
  colnames(obs_num) = c("Year", "Region", "Freq")
  wide_obs_num = dcast(obs_num, Region~Year, value = "Freq")  
## Using Freq as value column: use value.var to override.
  df = data.matrix(wide_obs_num[,2:45])   
  as.table(df)     
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## C    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    1    0    4
## F    0    0    0    2    3    2    0    0    0    0    0    1    0    1
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    3   10    2   13    4    4    3    5    2    5    2    3    5    3
## I    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## J    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## K    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## L    0    0    0    0    0    0    0    2    0    0    0    0    0    0
## M    0    0    0    0    0    0    0    0    0    0    0    1    0    0
## N    2    6    0    2    1    0    0    0    0    0    0    1    4    2
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## B    0    0    5    0    0    0    0    1    0    0    0    0    0    0
## C    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E   12    1    2    0    1    0    0   11    0    2    0    6    0    1
## F    0    1    1    0    1    0    0    0    2    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    5    2    3    1    1    0    4    1    4    4    0    1    1
## I    0    0    0    0    0    0    2    0    0    0    0    0   16   10
## J    0    0    1    1    1    0    3    1    1    0    0    0    2    2
## K    0    0    0    0    0    0    0    0    0    0    0    0    2    1
## L    0    0    0    0    0    0    0    0    0    1    1    4    4    2
## M    0    0    0    0    0    0    0    1    0    0    0    0    0    0
## N    0    1    5    1    2    0    1    1    1    1    1    1    1    1
## O    0    0    0    0    0    0    0    0    1    1    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    1    0    0    0    1    0    0    4    0    0    0
## C    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## D    0    0    0    0    0    0    0    0    0    0    0    4    0    1
## E    1    0    2    0    4    2    2    0    1    0    2    0    0    0
## F    0    0    1    1    0    1    3    6    2    2    0    0    1    0
## G    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## H    0    1    4   14   12    1    1    1    1    0    1    2    0    0
## I    0    0    0    6    3    0    1    1    4    0   17    6    2    0
## J    0    0    0    1    2    0    0    0    0    5    7    0    0    1
## K    1    5    5    3    2    0    0    0    0    0    2    0    0    0
## L    0    0    2    0    1    0    1    3    2    0   17    3    1    1
## M    0    0    0    0    0    1    2   14    8    0    0    0    0    1
## N    1    1    1    1    1    1    0    0    0    0    1    0    0    0
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    2    0    0    0    0    0    0    0    0
##   2013 2014
## A    0    0
## B    1    1
## C    2    1
## D    4    2
## E    0    0
## F    0    0
## G    0    0
## H    0    0
## I    2    2
## J    0    0
## K    0    0
## L    1    3
## M    0    0
## N    0    0
## O    0    0
## P    0    0
  df[ df>1 ] <- 1 
  levels(cover$Region) # How many regions? = how many colours
##  [1] "Clipperton"                 "Cocos_Islands"             
##  [3] "Colombia"                   "Costa_Rica_Central"        
##  [5] "Costa_Rica_Southern"        "Eastern_Galapagos_Islands" 
##  [7] "El_Salvador"                "Gulf_of_Chiriqui"          
##  [9] "Mexican_Central"            "Mexican_Northern"          
## [11] "Mexican_Southern"           "Nicaragua_Papagayo_Zone"   
## [13] "Northern_Galapagos_Islands" "Panama_Bight"              
## [15] "Revillagigedo_Islands"      "Western_Galapagos_Islands"
  Regions <- c(levels(obs_num$Region))
  color.regions = warmcool(17)
  
  # Figure 2a
  par(mfcol=c(2,1), cex.main = 1.5, mar = c(4,4,1,1), cex.lab = 1.5, cex.axis = 1,  las = 1)
  barplot(df,
          ylim=c(0,15),
          col = color.regions,   
          ylab="Number of regions surveyed per year",  
          xlab="Year",
          #lwd = 1,          
          las = 2,           
          axis.lty = 1, 
          cex.names = 0.4,   
          cex.axis = 0.8)
# Legend as a box
  legend( #"bottomright",    
     x = 2, y = 16, 
     xpd = NA,  # or TRUE 
     legend = Regions,          
     fill = color.regions[1:17],   
     inset=c(-0.5,0),
     cex = 0.35,
     text.font=0.5,            
     x.intersp = -0.5,
     y.intersp = 1,
     lwd = 0.1,
     bty="n")

# Figure 2b
obs_num = as.data.frame(table(cover$Year,cover$Country))    # Table per country
colnames(obs_num) = c("Year", "Country", "Freq")
wide_obs_num = dcast(obs_num, Country~Year, value = "Freq")   
## Using Freq as value column: use value.var to override.
dfc = data.matrix(wide_obs_num[,2:45])   
as.table(dfc)       
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## B    0    0    2    9    0    0    0    2    0    2    0    1    0    4
## C    0    0    0    2    3    2    0    0    0    0    0    2    0    1
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    5   16    0    6    5    4    3    5    2    4    2    4    9    5
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## B   12    1    7    0    1    0    0   15    0    6    4   10    4    3
## C    0    1    1    0    1    0    0    1    2    0    0    0    0    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## F    0    0    1    1    1    0    5    1    2    1    0    0   20   13
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    6    7    4    3    1    1    2    2    2    2    1    2    2
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## B    1    0    7    1    5    2    3    4    3    0   12    7    1    2
## C    0    0    1    1    0    4    5   20   10    2    0    0    1    1
## D    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    1    5    5   10    7    0    1    1    4    5   26    6    2    1
## G    0    0    0    0    0    0    0    0    0    0   11    0    0    0
## H    1    2    2   15   13    2    1    1    1    0    2    2    0    0
##   2013 2014
## A    2    1
## B    6    6
## C    0    0
## D    0    0
## E    0    0
## F    2    2
## G    0    0
## H    0    0
dfc[ dfc>1 ] <- 1 # replace all values greater than one with one
color.countries <- brewer.pal(9, "Oranges")
countries = c(levels(obs_num$Country))
barplot(dfc,
        col = color.countries,   
        ylim=c(0,10),
        ylab="Number of countries surveyed per year",  
        xlab="Year",
        las = 2,           
        axis.lty = 1,      
        cex.names = 0.4,   
        cex.axis = 0.8)
# Legend as a box
legend( #"bottomright",    
   x = 2, y = 10, 
   xpd = NA,  # or TRUE 
   legend = countries,         
   fill = color.countries[1:9],   
   inset=c(-0.5,0),
   cex = 0.35,
   text.font=1,            
   x.intersp = -0.5,
   y.intersp = 1,
   lwd = 0.1,
   bty="n")

FIGURE S3: Variation in the live coral cover for 1970–2014 per thermal regime

Temporal variation in the live coral cover for 1970–2014. * (a) Tropical upwelling regime. * (b) Seasonal regime. * (c) Thermally stable regime. * (d) Equatorial upwelling regime. The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.

# dev.off() # uncomment to reset parameters for plot

library(ggplot2)
library(lemon)

aggr.region$Thermal_regime<-factor(aggr.region$Thermal_regime, 
                                      levels=c("Upwelling", "Seasonal", "Thermally_Stable", "Galapagos_Islands"), 
                                      labels=c("(a)", "(b)", "(c)", "(d)"))
summary(aggr.region$Thermal_regime)
## (a) (b) (c) (d) 
##  75  28  74  25
ggplot(aggr.region, aes(x=Year, y=Coral_cover)) + 
  theme_bw() + theme(panel.border = element_blank(),   # set ggplot to white background
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(),
                              strip.background = element_blank(),
                              strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
  geom_point()+ geom_smooth (span = 0.3) +
  scale_x_continuous(limits=c(1970, 2014), expand = c(0.01, 0.01), "Year") +
  scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Coral cover") +
  facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top")